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It is frequently asserted that in a chaotic system two initially close points 
will separate at an exponential rate governed by the largest global Lyapunov 
exponent. Local Lyapunov exponents, however, are more directly relevant 
to predictability. The difference between the local and global Lyapunov ex- 
ponents, the large variations of local exponents over an attractor, and the 
saturation of error growth near the size of the attractor — all result in non- 
exponential scalings in errors at both short and long prediction times, some- 
times even obscuring evidence of exponential growth. Failure to observe ex- 
ponential error scaling cannot rule out deterministic chaos as an explanation. 
We demonstrate a simple model that quantitatively predicts observed error 
scaling from the local Lyapunov exponents, for both short and surprisingly 
long times. We comment on the relevance to atmospheric predictability as 
studied in the meteorological literature. 
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Chaotic behavior in dynamical systems is intermediate between precisely predictable, 
regular evolution and completely unpredictable, stochastic evolution. The global Lyapunov 
exponents [|1],|2| quantify the evolution of perturbations as the system evolves for long times. 
A positive global Lyapunov exponent implies that errors grow to the overall size of the 
attractor limiting predictability to finite times. 

Given observed data from a chaotic source, one can construct accurate empirical predic- 
tors Hf|J| which allow predictions for finite times — without any knowledge of the underlying 
dynamics. Positive global Lyapunov exponents causes the errors in these predictions to grow 
exponentially rapidly, and it is conventionally assumed that the prediction error E(t) will 
grow as 



Ai is the largest global Lyapunov exponent. Indeed this error growth is used as a diagnostic 
of deterministic chaos in the analysis of observed data j|,[5] . 

This paper is devoted to a closer examination of this common assumption, and we will 
show it is often incorrect for times where the system is predictable. For t — > oo the largest 
global Lyapunov exponent Ai is defined in terms of a long-term average growth rate, but 
it does not always reflect the actual growth of prediction errors in finite time. Further for 
t — > oo the perturbed orbit, although eventually uncorrelated with the reference orbit, is 
constrained to stay on the attractor, and thus the size of a perturbation will saturate near 
the overall size of the attractor. Reconciling these two aspects of error growth leads to our 
improvement to Equation ([!]). 

We briefly review the definition of local Lyapunov exponents. Our discrete time dynam- 
ical system is y(n + 1) = F(y(n)). Small perturbations 5y(n) to this orbit evolve L steps 
forward in time according to the linearized dynamics 



With DF a fc(y) = dF a (y) /dyb the Jacobian matrix of the dynamics, we denote DF L (y(n) = 
Uf=o DF(y(n + i))- The Oseledec matrix §, 



E(t) = £(0)exp(Ait). 



(1) 




(2) 




(3) 
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has eigenvalues e Xl( ^ ,L \ e X2< ^ ,L \ . . . , e Xd ^ ,L \ in d-dimensions. We order the exponents as 
Ai > A 2 > A 3 . . . > Xd- These local exponents A a (y, L) address the growth (or decay) over 
L time steps of perturbations made around some point y in the state space. As L — > 00 
X a (y,L) — > A a , the global exponents. 

Taken over an initially isotropic distribution the initial root-mean-squared Euclidean 
distance of the evolved perturbation vectors to the origin will grow by Yh=\ A 2 ) 1 / 2 - 
A a (x, L) = exp[LA a (x, L)}. Lorenz first derived this formula in an early paper ||. We 
employ the analogous quantity appropriate for a geometric mean instead of the arithmetic 
mean. We found this to be \ Sf=i A* via numerical experiment, though we do not yet know 
of a simple, general derivation corresponding to that in reference || . We note that direct nu- 
merical computation of the matrix product of many Jacobians leads to ill-conditioning, and 
thus for practical computation we employ the algorithm of reference || to stably compute 
the local Lyapunov exponents. 

The quantity E(x,L) = ~ Yh=i Aj(x, L), which we denote the "expansion factor", quan- 
tifies the multiplicative growth of typical predictor errors starting at x, looking ahead L 
steps. The expansion factor is best calculated from local exponents as 

log E(x, L) = LAx(x, L) + log[l + e ^(x,L)-A l( x,L)) + ^(x^-mx.l)) + 

+ e L(A d (X,L)-A 1 (X,L)) ] _ logrf (4) 

For increasing L, Ai(x, L) quickly dominates the expansion factor. This also motivates Equa- 
tion ([[]), but here Ai(x, L), is the finite-time Lyapunov exponent, which does not converge 
very quickly to Ai. 0. 

The other ingredient in our model for the average prediction error is a saturation cutoff. 
If we limit the maximal expansion factor for each initial condition to a constant R, we then 
compute the geometric mean (over reference points on the attractor x(i)) of the expansion 
factors -E"(x, L), which is the arithmetic mean of log-E^x, L), hard-limited by p = log-R: 

1 N 

\ogX(L,p) = -]Tmm[log£(x(*),L),p] -p. (5) 

iV i=l 

The saturation cutoff models the fact that finite-sized perturbations and prediction errors 
cannot grow indefinitely. We estimate R is the ratio of the geometric mean of |x(i) — x(j)| 
over uncorrelated pairs of attractor points to the geometric mean of the initial perturbation 
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magnitude. Larger R corresponds to better predictability. Note that we are averaging 
individually thresholded expansion factors — not thresholding an average expansion factor. 
Before the saturation sets in \ogX(L) ~ LX(L), with X(L) the average finite-time Lyapunov 
exponent [[F]|§. 

We now compare (1) the growth of errors actually incurred by making repeated predic- 
tions and (2) the growth rate implied by the thresholded expansion factors, Equation (^|). 
This suggests that prediction errors grow as do the separation of initially close trajectories. 
This means we are considering the errors in prediction that arise from inherent dynamical 
instability and not inaccuracies as a result of specific features of the prediction scheme. 

Our measure of average prediction error is the geometric mean over initial conditions of an 
L-step iterated predictor, a composition of L one-step predictors. The error is normalized by 
the geometric mean distance between all time-decorrelated pairs of points on the attractor, 
so that the absence of predictability corresponds to zero logarithmic error. With G(«) the 
one-step predictor, the normalized prediction error L steps ahead is 

i N i N 

logx(L) = - £log |G L (x(z)) - x(i + L)\ - -j £ log |x(i) - x(j)|. (6) 
iV t=i iy i,j=i 

As for \ogX(L, p) liniL-KX) \ogx{L) = 0. The iterated prediction is not rescaled to remain 
close to the reference trajectory. 

The main empirical result is that on chaotic attractors \ogx{L) = logX(L, p) given the 
correct threshold p. We evaluate the local Lyapunov exponents, then compute the single 
parameter family of curves given by Equation (|5|) over a range of p. We find that once 
the best value of p is found, the resulting logX(L,p) quantitatively predicts the scaling 
of errors throughout the range of time examined. That equation (|5|) predicts the scaling of 
errors in the saturation region (L — > oo) is somewhat surprising because Lyapunov exponents 
quantify linear expansion rates of infinitesimal perturbations, but in the saturation region 
one envisions typically large deviations. A possible explanation is that in the saturation 
regime average prediction error is controlled by the tail of the distribution of individual 
expansion factors: those initial conditions that happen to be exceptionally predictable and 
whose expansion factors thus remain below the threshold for especially long times. This is 
plausible because local Lyapunov exponents often have wide distributions 

Figure |T] shows iterated prediction errors,log 10 (x(£)), computed using nonlinear kernel 
regression || versus thresholded expansion factors (|]) for data from the Ikeda map of the 
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plane to itself fL0| . The data set is 20,000 real and imaginary components from the complex- 
valued map z(n + 1) = p + Bz(n) exp[in — ia/(l + |z(^)| 2 )] with parameters p = 1.0, B = 
0.9, k = 0.4, a = 6.0. At short times, there is an obvious exponential scaling region (a 
line with a constant slope Ai), as predicted by Equation ([[J). Starting at L ~ 10 the slope 
decreases and curves off to zero, well modeled by the thresholded expansion factor. We did 
not optimize p in any sophisticated manner, but simply stepped p by 0.2 and selected the best 
match. From a time series y(i),i = 1 ... N, we compute a prediction for an input vector x as 
G(x) = (Eili y(* + l)#(x - y(i))) I (Eili K(x - y(i))) using a Gaussian kernel K(z) = 
exp(— \z\ 2 /a 2 ) with o set to twice the geometric mean nearest-neighbor distance of the y(z). 
This is a global prediction formula but effectively functions as a localized interpolator. The 
local Lyapunov exponents were calculated from this same data set || without the equations. 

Figure |2| shows the same information for data from a three dimensional flow of Lorenz [JTT . 
The dynamical equations are x = —y 2 — z 2 — a(x — F),y = xy — bxz — y + G,z = bxy + xz — z, 
with parameters set a = 0.25, b = 4.0, F = 8.0, and G = 1.0, sampled at intervals St = 0.2. 
The attractor has a dimension of about 2.5. There is no obvious exponential scaling regime 
here. Given only the curve of predictor errors, one could not easily identify the Lyapunov 
exponent, in contrast to the previous example. The appropriately thresholded expansion 
factor, however, agrees with the observed scaling of prediction errors. For contrast, the 
graph also shows the results of choosing either too small a threshold (the circles) or too 
large a threshold (the triangles) for this prediction scheme and data set. 

Now we examine our main result logx(L) = log X(L,p) with a better controlled ex- 
periment. This time as our "prediction function" we use the actual dynamical equations 
of the Lorenz system but start the integration with an initial condition slightly perturbed 
from the reference point by a small uniformly distributed random vector t](i): G L (x(i)) = 
F L (x(i) +Tf(i)). We directly calculate Lyapunov exponents from the known differential equa- 
tions by simultaneously integrating the equations of motion and the tangent space flow. Fig- 
ure ^| compares logx(^) with \ogX(L, p) with varying sizes of the initial perturbation. The 
growth of the perturbations matches the thresholded expansion factor, with the previously 
free parameter p fixed at the predicted value p = (log |x(z) — X-(J)\)i,j=i„jf— (log \7}(i)\)i=i,„N> 
confirming our model. 

This scaling with L is generic to most low dimensional flows: a curve at low L scaling with 
an initially high slope, because the average largest local exponent is larger than the global 
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exponent ||, an exponential region (constant slope) corresponding to the global Lyapunov 
exponent, and a long tail curving off to zero as the finite-size threshold takes effect. The 
size and existence of the region of constant slope depends on the size of the initial error. If 
the error is large enough, there may be no region of exponential scaling as the curvature due 
to finite-time exponents merges with the curvature due to thresholds. Higher-dimensional 
data are likely to have lower initial predictability than very low- dimensional attractors, and 
therefore likely to show little observably exponential error scaling. Intermittent dynamics 
will create a wide variation in finite-time Lyapunov exponents, thus moving forward the 
time where saturation begins to take effect. In analyzing dynamical systems more complex 
than simple one or two dimensional models, requiring manifestly exponential error scaling 
as confirming deterministic chaos is unrealistic. 

We have identified empirical prediction error with the evolution of perturbations, but 
the identity does not always strictly hold. Figure ^ shows the prediction errors on the same 
Lorenz flow |Tl| data using an accurate local quadratic polynomial technique H,§|. This 
time, the prediction error is not satisfactorily matched by the thresholded expansion factor 
with any p, the main difference being a substantially larger slope at small L. We attribute 
part of this discrepancy to the assumption, used in deriving the expansion factor, that the 
initial perturbations are distributed uniformly in all directions. When a very accurate model, 
such as this one, is combined with extremely clean data, short term forecasting errors will 
be quite small, and the predicted trajectory will closely shadow the attractor, and thus, the 
unstable manifold. Therefore the initial rate of expansion will be larger than for isotropically 
distributed errors, being described better by a new "primary" expansion factor LAi(x, L) 
that only considers error expansion due to the single largest local Lyapunov exponent. This 
quantity increases faster at short times than the standard expansion factor, but still, it fails 
to precisely match the scaling of forecast error. Another discrepancy not accounted for by 
the expansion factor is that iterating imperfect models injects new error at every timestep, 
not only at the beginning. Usually exponential expansion of old error dominates new error, 
except at the shortest times. This effect will cause yet another increase to the slope at 
small times. When we add some isotropic noise to the initial condition before applying the 
iterated local predictor, outstanding agreement with the scaling predicted by Equation ([5]) 
is restored. We note that the divergence between straight prediction error and expansion 
factor is exaggerated by the particular conditions in force here: very good prediction on a 
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data base of very clean low-dimensional data from a smooth chaotic flow. We performed 
the same test on noisier experimental chaotic data and saw a smaller difference. 

Is adding noise to the initial condition "cheating" ? To be completely rigorous, one ought 
to model the distribution of prediction errors in the various directions corresponding to the 
local Lyapunov exponents. This is not possible without knowledge of detailed properties 
of the specific prediction scheme, and is obviously beyond the scope of this paper. The 
approximation of isotropic perturbation may often be acceptable, and can be enforced by 
adding artificial noise. In the context of analyzing observed data, a successful match be- 
tween prediction errors (even with artificial initial perturbations) and thresholded expansion 
factors is a novel cross-check that affirms the validity of the modeling procedures, even if 
one cannot exactly quantify the error scaling for a particular prediction scheme. Accurate 
evaluation of Lyapunov exponents from degraded or high- dimensional data is somewhat 
difficult in practice, requiring good estimates of derivatives of the implied evolution func- 
tion. Various reconstruction parameters, such as embedding and local manifold dimensions 
and time delay |8||| may yield substantially different numerical answers for Lyapunov ex- 
ponents, even if all are topologically acceptable. Further requiring a good match between 
expansion factors and error scaling provides a criterion for selecting among the otherwise 
equivalent choices. In general, we have found that prediction error scaling varies less with 
reconstruction parameters than Lyapunov exponents. 

Figure |5| shows (log£/(x(z), L)) (logX without the threshold) for the Lorenz flow. At 
the smallest times, the expansion factor does not in fact possess a constant slope, and only 
flattens out to a straight line, with slope equal to the infinite-time Lyapunov exponent after 
an initial transient regime. This curvature is always in the direction seen on this graph 
(higher slope at shorter times) and is a result of the fact that the average local Lyapunov 
exponent X(L) approaches the global exponent from above as L — > oo ||. This curvature 
explains non-exponential scaling of prediction error for early time intervals. Also shown is 
the mean plus and minus one standard deviation of the distribution of individual expansion 
factors E(pc, L). The Figure demonstrates that the variation in local exponents as a function 
of initial condition causes a wide variation in expansion factor that increases with increasing 
L. The interaction of this wide distribution with the threshold results in the saturation 
of prediction error for longer times. The effects of the finite size of the attractor begin 
to manifest themselves when an appreciable number of the individual expansion factors 
approach the cutoff, which occurs well before the mean expansion factor does so. 
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Notice that the width of the distribution increases as L — > oo; we have observed this 
feature in all systems we examined. An important conclusion is that even though the 
local Lyapunov exponent converges to a single global exponent which is independent of 
initial condition, i.e. lim^oo Ai(x, L) = Ai, the expansion factors of various individual 
initial conditions do not do so: lim^oo E(pc, L) ^ E(L). This is because the standard 
deviation of the distribution of local exponents ((A(x, L) — A(L)) 2 ) typically decreases at 
a rate substantially slower than L~ x . Of course, considering global saturation, normalized 
error will eventually converge to unity, but this is an entirely different mechanism. The 
message is that considering dynamical predictability solely in view of the largest global 
Lyapunov exponent may be conceptually misleading as well as quantitatively inaccurate. 

The temporal development of forecast error has been a prime concern of the weather 
forecasting community for many years, starting with work of Lorenz ||12|. Common practice 



has been to hypothesize ad hoc relationships for the average error as a function of time; 



usually, in the form of a differential equation. Lorenz found |T2[ that increase of mean 
deviation between initially close initial atmospheric states could be fitted by an empirical 
law of the form 

E(t) = aE(t)(l-E(t)/E(oo)). (7) 

A positive Lyapunov exponent motivates the first term; the inevitable saturation at max- 
imum error, the second. Stroe and Royer [13[] compared generalized parameterizations 



of empirical growth laws that include Lorenz's with results of large-scale atmospheric 
simulations and some experimental observations. The scaling of mean error at larger 
times, in the saturation regime, appeared to be better fit by an exponential law such as 
dlogE(t)/dt = — (3 log E(t), similar to results seen in our work, but a single empirical rule 
governing the initial growth of errors was not clearly indicated, either in this or previ- 
ous studies. We explain this with the fact that the initial growth of error is governed by 
finite-time, rather than global, Lyapunov exponents. This results in an initial regime of 
non-exponential growth, second, this error growth rate depends on the coordinate system 
used 0. The differing measures of phase-space distance employed by atmospheric scientists 
will result in different growth rates making the notion of a universal "doubling time for small 
errors" less useful than commonly believed. Our modeling of the error growth, which holds 
the mean finite-time Lyapunov exponents responsible at small error, but with their variation 
most important at saturation, backs up Stroe and Royer's conclusions that the infinitesimal 
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growth rate cannot be easily deduced from the saturation rate via fitting a single empirical 
formula like equation (|7]) . One further point that we wish to make is that we have observed 
that using the arithmetic mean for the ensemble average of both error and expansion fac- 
tors rather than the geometric results in far more "noisy" curves, and thus is it not clear 
whether expansion factors accurately match prediction errors with ensemble averaging in the 
arithmetic sense, though we suspect so. The arithmetic mean, commonly employed in me- 
teorological literature, seems to be dominated by fluctuations in a few samples in the large 
error tail of the distribution. Convergence of the arithmetical ensemble average appears 
to require numbers of points excessively large even for this low-dimensional investigation. 
Choosing an arithmetical average also appears to accentuate the initial superexponential 
growth. Still, the qualitative behavior of error scaling seen in large-scale atmospheric simu- 
lations and experiment |13| appears compatible with our model, which we suggest as a more 
fundamental explanation for observed error growth. In the nonlinear dynamics literature, 
one example in Farmer and Sidorowich Q demonstrated initially faster-than-exponential 
error scaling, but remained unexplained in that work. 

With the exception of the work of Lorenz [[J , the direct use of local Lyapunov exponents 
to quantify error growth in atmospheric dynamics is rather recent fll4|-|l6f and remains open 
to further development. The Lyapunov exponents have generally been previously considered 
only in the context of infinitesimal errors. This present work shows that accounting for a 
threshold apparently allows finite-time Lyapunov exponents to quantify error scaling at 
both small and substantially larger levels of error. Work remains concerning more realistic 
models, of course. Results of other low- dimensional chaotic data sets that we have examined, 
including experimental data, agree with the conclusions of this paper. 
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FIGURES 



FIG. 1. log 10 x(X) (solid line) and log 10 X(L,p) (circles) for the Ikeda map. p = 2.4 gives a 
good fit. 



FIG. 2. log 10 x(L) (solid line) and log 10 X(L,p) for data from the Lorenz flow, with p = 1.2 
(triangles), p = 1.6 (circles), and p = 2.0 (diamonds). 



FIG. 3. log 10 x(-k) (solid lines) and log 10 X(L,p) (symbols) using known flow equations and 
varying sizes of initial perturbations. 

FIG. 4. log 10 x(-^) and log 10 X(L,p) (symbols) using local quadratic prediction, with no initial 
perturbation (solid line), and a 1% perturbation (dashed line) before prediction. 

FIG. 5. Geometric mean of local expansion factors (solid line), plus one standard deviation of 
the distribution (dotted line) and minus one standard deviation (dashed line), for Lorenz flow. 
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Chaotic behavior in dynamical systems is intermediate between precisely predictable, 
regular evolution and completely unpredictable, stochastic evolution. The global Lyapunov 
exponents [1,2] quantify the evolution of perturbations as the system evolves for long times. 
A positive global Lyapunov exponent implies that errors grow to the overall size of the 
attractor limiting predictability to finite times. 

Given observed data from a chaotic source, one can construct accurate empirical predic- 
tors [3,4,6] which allow predictions for finite times — without any knowledge of the underlying 
dynamics. Positive global Lyapunov exponents causes the errors in these predictions to grow 
exponentially rapidly, and it is conventionally assumed that the prediction error Eit) will 
grow as 



Ai is the largest global Lyapunov exponent. Indeed this error growth is used as a diagnostic 
of deterministic chaos in the analysis of observed data [4,5]. 

This paper is devoted to a closer examination of this common assumption, and we will 
show it is often incorrect for times where the system is predictable. For t — > oo the largest 
global Lyapunov exponent Ai is defined in terms of a long-term average growth rate, but 
it does not always reflect the actual growth of prediction errors in finite time. Further for 
t — > oo the perturbed orbit, although eventually uncorrelated with the reference orbit, is 
constrained to stay on the attractor, and thus the size of a perturbation will saturate near 
the overall size of the attractor. Reconciling these two aspects of error growth leads to our 
improvement to Equation (1). 

We briefly review the definition of local Lyapunov exponents. Our discrete time dynam- 
ical system is y(n + 1) = F(y(n)). Small perturbations 6y(n) to this orbit evolve L steps 
forward in time according to the linearized dynamics 







With DF a fe(y) = dF a (y) / dy\, the Jacobian matrix of the dynamics, we denote DF L (y(n) 
Ylfjo DF(y(n + i)). The Oseledec matrix [2], 




(3) 
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has eigenvalues e Xl ^' L \ e^ 2 ^ ,L \ . . . , e Xd ^' L \ in d- dimensions. We order the exponents as 
Ai > A 2 > A 3 . . . > Arf. These local exponents A a (y, L) address the growth (or decay) over 
L time steps of perturbations made around some point y in the state space. As L — > oo 
X a (y } L) — > A a , the global exponents. 

Taken over an initially isotropic distribution the initial root-mean-squared Euclidean 
distance of the evolved perturbation vectors to the origin will grow by J2i=i A. 2 ) 1 / 2 . 
A a (x, L) = exp[_LA a (x, L)]. Lorenz first derived this formula in an early paper [9]. We 
employ the analogous quantity appropriate for a geometric mean instead of the arithmetic 
mean. We found this to be ^ J2i=i via numerical experiment, though we do not yet know 
of a simple, general derivation corresponding to that in reference [9]. We note that direct nu- 
merical computation of the matrix product of many Jacobians leads to ill-conditioning, and 
thus for practical computation we employ the algorithm of reference [8] to stably compute 
the local Lyapunov exponents. 

The quantity E(x, L) = ^ J2i=\ A 4 -(x, L), which we denote the "expansion factor", quan- 
tifies the multiplicative growth of typical predictor errors starting at x, looking ahead L 
steps. The expansion factor is best calculated from local exponents as 

log£(x, L) = LX 1 (^,L) + log[l + e £(MX,L)-A l( X,L)) + e L(A 3 (X,L)-A l( X,L)) + _ _ _ 

+ e L(A d (X,L)-A 1 (X,L)) ] _ logJ (4) 

For increasing L } Ai(x, L) quickly dominates the expansion factor. This also motivates Equa- 
tion (1), but here Ai(x,iv), is the finite-time Lyapunov exponent, which does not converge 
very quickly to Ai. [7]. 

The other ingredient in our model for the average prediction error is a saturation cutoff. 
If we limit the maximal expansion factor for each initial condition to a constant i?, we then 
compute the geometric mean (over reference points on the attractor x(z)) of the expansion 
factors -E(x, L), which is the arithmetic mean of log -E(x, L), hard-limited by p = log R: 

1 N 

log X(L, p) = — min [log £(x(z), L),p\- p. (5) 

8 = 1 

The saturation cutoff models the fact that finite-sized perturbations and prediction errors 
cannot grow indefinitely. We estimate R is the ratio of the geometric mean of |x(z) — x(j)| 
over uncorrelated pairs of attractor points to the geometric mean of the initial perturbation 
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magnitude. Larger R corresponds to better predictability. Note that we are averaging 
individually thresholded expansion factors — not thresholding an average expansion factor. 
Before the saturation sets in logX(L) LX(L) } with X(L) the average finite-time Lyapunov 
exponent [7,8]. 

We now compare (1) the growth of errors actually incurred by making repeated predic- 
tions and (2) the growth rate implied by the thresholded expansion factors, Equation (5). 
This suggests that prediction errors grow as do the separation of initially close trajectories. 
This means we are considering the errors in prediction that arise from inherent dynamical 
instability and not inaccuracies as a result of specific features of the prediction scheme. 

Our measure of average prediction error is the geometric mean over initial conditions of an 
Z-step iterated predictor, a composition of L one-step predictors. The error is normalized by 
the geometric mean distance between all time-decorrelated pairs of points on the attractor, 
so that the absence of predictability corresponds to zero logarithmic error. With G(«) the 
one-step predictor, the normalized prediction error L steps ahead is 

i N -r N 

lo S X(L) = - £log |G L (x(z)) - x(i + L)\- — £ log |x(i) - x(j)|. (6) 

i=l 2J = 1 

As for log X(L, p) lim^^oo log = 0. The iterated prediction is not rescaled to remain 

close to the reference trajectory. 

The main empirical result is that on chaotic attractors logx(L) = log X(L, p) given the 
correct threshold p. We evaluate the local Lyapunov exponents, then compute the single 
parameter family of curves given by Equation (5) over a range of p. We find that once 
the best value of p is found, the resulting log X(L } p) quantitatively predicts the scaling 
of errors throughout the range of time examined. That equation (5) predicts the scaling of 
errors in the saturation region (L — > oo) is somewhat surprising because Lyapunov exponents 
quantify linear expansion rates of infinitesimal perturbations, but in the saturation region 
one envisions typically large deviations. A possible explanation is that in the saturation 
regime average prediction error is controlled by the tail of the distribution of individual 
expansion factors: those initial conditions that happen to be exceptionally predictable and 
whose expansion factors thus remain below the threshold for especially long times. This is 
plausible because local Lyapunov exponents often have wide distributions [7,8]. 

Figure 1 shows iterated prediction errors, log 10 (x(iv)), computed using nonlinear kernel 
regression [6] versus thresholded expansion factors (5) for data from the Ikeda map of the 
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plane to itself [10]. The data set is 20,000 real and imaginary components from the complex- 
valued map z(n + 1) = p + Bz(n) exp[in — ia/(l + \z(n) | 2 )] with parameters p = 1.0, B = 
0.9, k = 0.4, a = 6.0. At short times, there is an obvious exponential scaling region (a 
line with a constant slope Ai), as predicted by Equation (1). Starting at L 10 the slope 
decreases and curves off to zero, well modeled by the thresholded expansion factor. We did 
not optimize p in any sophisticated manner, but simply stepped p by 0.2 and selected the best 
match. From a time series y(i), i = 1 . . . N, we compute a prediction for an input vector x as 
G(x) = (j2iLi y(i + l)A"(x — y(z))) / (j2iLi A"( x _ y(0)) usm g a Gaussian kernel K(z) = 
exp( — |z| 2 /cr 2 ) with a set to twice the geometric mean nearest-neighbor distance of the y(i). 
This is a global prediction formula but effectively functions as a localized interpolator. The 
local Lyapunov exponents were calculated from this same data set [8] without the equations. 

Figure 2 shows the same information for data from a three dimensional flow of Lorenz [11]. 
The dynamical equations are x = —y 2 — z 2 — a(x — F) } y = xy — bxz — y-\-G } z= bxy -\-xz — z } 
with parameters set a = 0.25, b = 4.0, F = 8.0, and G = 1.0, sampled at intervals St = 0.2. 
The attractor has a dimension of about 2.5. There is no obvious exponential scaling regime 
here. Given only the curve of predictor errors, one could not easily identify the Lyapunov 
exponent, in contrast to the previous example. The appropriately thresholded expansion 
factor, however, agrees with the observed scaling of prediction errors. For contrast, the 
graph also shows the results of choosing either too small a threshold (the circles) or too 
large a threshold (the triangles) for this prediction scheme and data set. 

Now we examine our main result logx(iv) = log X(L } p) with a better controlled ex- 
periment. This time as our "prediction function" we use the actual dynamical equations 
of the Lorenz system but start the integration with an initial condition slightly perturbed 
from the reference point by a small uniformly distributed random vector G L (x(z)) = 

F L (x(z) -\-t](i)). We directly calculate Lyapunov exponents from the known differential equa- 
tions by simultaneously integrating the equations of motion and the tangent space flow. Fig- 
ure 3 compares log x{L) with log X(L } p) with varying sizes of the initial perturbation. The 
growth of the perturbations matches the thresholded expansion factor, with the previously 
free parameter p fixed at the predicted value p = (log |x(z) — x(j) |)i J= i...jv — (log \ f]{i)\)i=\...Ni 
confirming our model. 

This scaling with L is generic to most low dimensional flows: a curve at low L scaling with 
an initially high slope, because the average largest local exponent is larger than the global 
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exponent [8], an exponential region (constant slope) corresponding to the global Lyapunov 
exponent, and a long tail curving off to zero as the finite-size threshold takes effect. The 
size and existence of the region of constant slope depends on the size of the initial error. If 
the error is large enough, there may be no region of exponential scaling as the curvature due 
to finite-time exponents merges with the curvature due to thresholds. Higher-dimensional 
data are likely to have lower initial predictability than very low-dimensional attractors, and 
therefore likely to show little observably exponential error scaling. Intermittent dynamics 
will create a wide variation in finite-time Lyapunov exponents, thus moving forward the 
time where saturation begins to take effect. In analyzing dynamical systems more complex 
than simple one or two dimensional models, requiring manifestly exponential error scaling 
as confirming deterministic chaos is unrealistic. 

We have identified empirical prediction error with the evolution of perturbations, but 
the identity does not always strictly hold. Figure 4 shows the prediction errors on the same 
Lorenz flow [11] data using an accurate local quadratic polynomial technique [4,5]. This 
time, the prediction error is not satisfactorily matched by the thresholded expansion factor 
with any p } the main difference being a substantially larger slope at small L. We attribute 
part of this discrepancy to the assumption, used in deriving the expansion factor, that the 
initial perturbations are distributed uniformly in all directions. When a very accurate model, 
such as this one, is combined with extremely clean data, short term forecasting errors will 
be quite small, and the predicted trajectory will closely shadow the attractor, and thus, the 
unstable manifold. Therefore the initial rate of expansion will be larger than for isotropically 
distributed errors, being described better by a new "primary" expansion factor L\i(x.,L) 
that only considers error expansion due to the single largest local Lyapunov exponent. This 
quantity increases faster at short times than the standard expansion factor, but still, it fails 
to precisely match the scaling of forecast error. Another discrepancy not accounted for by 
the expansion factor is that iterating imperfect models injects new error at every timestep, 
not only at the beginning. Usually exponential expansion of old error dominates new error, 
except at the shortest times. This effect will cause yet another increase to the slope at 
small times. When we add some isotropic noise to the initial condition before applying the 
iterated local predictor, outstanding agreement with the scaling predicted by Equation (5) 
is restored. We note that the divergence between straight prediction error and expansion 
factor is exaggerated by the particular conditions in force here: very good prediction on a 
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data base of very clean low- dimensional data from a smooth chaotic flow. We performed 
the same test on noisier experimental chaotic data and saw a smaller difference. 

Is adding noise to the initial condition "cheating"? To be completely rigorous, one ought 
to model the distribution of prediction errors in the various directions corresponding to the 
local Lyapunov exponents. This is not possible without knowledge of detailed properties 
of the specific prediction scheme, and is obviously beyond the scope of this paper. The 
approximation of isotropic perturbation may often be acceptable, and can be enforced by 
adding artificial noise. In the context of analyzing observed data, a successful match be- 
tween prediction errors (even with artificial initial perturbations) and thresholded expansion 
factors is a novel cross-check that affirms the validity of the modeling procedures, even if 
one cannot exactly quantify the error scaling for a particular prediction scheme. Accurate 
evaluation of Lyapunov exponents from degraded or high-dimensional data is somewhat 
difficult in practice, requiring good estimates of derivatives of the implied evolution func- 
tion. Various reconstruction parameters, such as embedding and local manifold dimensions 
and time delay [8,3] may yield substantially different numerical answers for Lyapunov ex- 
ponents, even if all are topologically acceptable. Further requiring a good match between 
expansion factors and error scaling provides a criterion for selecting among the otherwise 
equivalent choices. In general, we have found that prediction error scaling varies less with 
reconstruction parameters than Lyapunov exponents. 

Figure 5 shows (log _E(x(z), L)) (logX without the threshold) for the Lorenz flow. At 
the smallest times, the expansion factor does not in fact possess a constant slope, and only 
flattens out to a straight line, with slope equal to the infinite-time Lyapunov exponent after 
an initial transient regime. This curvature is always in the direction seen on this graph 
(higher slope at shorter times) and is a result of the fact that the average local Lyapunov 
exponent X(L) approaches the global exponent from above as L — > oo [8]. This curvature 
explains non-exponential scaling of prediction error for early time intervals. Also shown is 
the mean plus and minus one standard deviation of the distribution of individual expansion 
factors -E(x, L). The Figure demonstrates that the variation in local exponents as a function 
of initial condition causes a wide variation in expansion factor that increases with increasing 
L. The interaction of this wide distribution with the threshold results in the saturation 
of prediction error for longer times. The effects of the finite size of the attractor begin 
to manifest themselves when an appreciable number of the individual expansion factors 
approach the cutoff, which occurs well before the mean expansion factor does so. 
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Notice that the width of the distribution increases as L — > oo; we have observed this 
feature in all systems we examined. An important conclusion is that even though the 
local Lyapunov exponent converges to a single global exponent which is independent of 
initial condition, i.e. lim^oo Ai(x, L) = Ai, the expansion factors of various individual 
initial conditions do not do so: lim^oo E(x., L) ^ E(L). This is because the standard 
deviation of the distribution of local exponents ((A(x,Z) — X(L)) 2 ) typically decreases at 
a rate substantially slower than L~ x . Of course, considering global saturation, normalized 
error will eventually converge to unity, but this is an entirely different mechanism. The 
message is that considering dynamical predictability solely in view of the largest global 
Lyapunov exponent may be conceptually misleading as well as quantitatively inaccurate. 

The temporal development of forecast error has been a prime concern of the weather 
forecasting community for many years, starting with work of Lorenz [9,12] . Common practice 
has been to hypothesize ad hoc relationships for the average error as a function of time; 
usually, in the form of a differential equation. Lorenz found [12] that increase of mean 
deviation between initially close initial atmospheric states could be fitted by an empirical 
law of the form 



A positive Lyapunov exponent motivates the first term; the inevitable saturation at max- 
imum error, the second. Stroe and Royer [13] compared generalized parameterizations 
of empirical growth laws that include Lorenz's with results of large-scale atmospheric 
simulations and some experimental observations. The scaling of mean error at larger 
times, in the saturation regime, appeared to be better fit by an exponential law such as 
dlog E(t) / dt = — f3 log E(t), similar to results seen in our work, but a single empirical rule 
governing the initial growth of errors was not clearly indicated, either in this or previ- 
ous studies. We explain this with the fact that the initial growth of error is governed by 
finite-time, rather than global, Lyapunov exponents. This results in an initial regime of 
non-exponential growth, second, this error growth rate depends on the coordinate system 
used [7]. The differing measures of phase-space distance employed by atmospheric scientists 
will result in different growth rates making the notion of a universal "doubling time for small 
errors" less useful than commonly believed. Our modeling of the error growth, which holds 
the mean finite-time Lyapunov exponents responsible at small error, but with their variation 
most important at saturation, backs up Stroe and Royer's conclusions that the infinitesimal 




(7) 
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growth rate cannot be easily deduced from the saturation rate via fitting a single empirical 
formula like equation (7). One further point that we wish to make is that we have observed 
that using the arithmetic mean for the ensemble average of both error and expansion fac- 
tors rather than the geometric results in far more "noisy" curves, and thus is it not clear 
whether expansion factors accurately match prediction errors with ensemble averaging in the 
arithmetic sense, though we suspect so. The arithmetic mean, commonly employed in me- 
teorological literature, seems to be dominated by fluctuations in a few samples in the large 
error tail of the distribution. Convergence of the arithmetical ensemble average appears 
to require numbers of points excessively large even for this low-dimensional investigation. 
Choosing an arithmetical average also appears to accentuate the initial superexponential 
growth. Still, the qualitative behavior of error scaling seen in large-scale atmospheric simu- 
lations and experiment [13] appears compatible with our model, which we suggest as a more 
fundamental explanation for observed error growth. In the nonlinear dynamics literature, 
one example in Farmer and Sidorowich [4] demonstrated initially faster-than-exponential 
error scaling, but remained unexplained in that work. 

With the exception of the work of Lorenz [9], the direct use of local Lyapunov exponents 
to quantify error growth in atmospheric dynamics is rather recent [14-16] and remains open 
to further development. The Lyapunov exponents have generally been previously considered 
only in the context of infinitesimal errors. This present work shows that accounting for a 
threshold apparently allows finite-time Lyapunov exponents to quantify error scaling at 
both small and substantially larger levels of error. Work remains concerning more realistic 
models, of course. Results of other low-dimensional chaotic data sets that we have examined, 
including experimental data, agree with the conclusions of this paper. 
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FIGURES 

FIG. 1. log 10 x(i) (solid line) and log 10 X(L, p) (circles) for the Ikeda map. p = 2.4 gives a 
good fit. 

FIG. 2. log 10 x(i) (solid line) and \og w X(L,p) for data from the Forenz flow, with p = 1.2 
(triangles), p = 1.6 (circles), and p = 2.0 (diamonds). 

FIG. 3. log 10 x(i) (solid lines) and log 10 X(L, p) (symbols) using known flow equations and 
varying sizes of initial perturbations. 

FIG. 4. log 10 x(i) and \og w X(L,p) (symbols) using local quadratic prediction, with no initial 
perturbation (solid line), and a 1% perturbation (dashed line) before prediction. 

FIG. 5. Geometric mean of local expansion factors (solid line), plus one standard deviation of 
the distribution (dotted line) and minus one standard deviation (dashed line), for Forenz flow. 
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